Diffusion Mechanisms in Lithium Disihcate Melt by Molecular Dynamics Simulation 
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In this work we study the difTusion mechanisms in lithium disilicate melt using molecular dynamics 
simulation, which has an edge over other simulation methods because it can track down actual atomic 
rearrangements in materials once a realistic interaction potential is applied. Our simulation results 
of diffusion coefficients show an excellent agreement with experiments. We also demonstrate that 
our system obeys the famous Stokes-Einstein relation at least down to 1400 K, while a decoupling 
between relaxation and viscosity takes place at a higher temperature. Additionally, an analysis 
on the dynamical behavior of slow-diffusing atoms reveals explicitly the presence of dynamical 
heterogeneities. 



I. INTRODUCTION 

The dynamics of complex systems has challenged sci- 
entists for many years since the first attempts on explain- 
ing the glass transition^"*. Glasses are among the most 
important materials of technological interest and thor- 
ough understanding of their synthesis is crucial. It is a 
well-known fact that diffusion controls processes like nu- 
cleation, growth, and electrical conductivity in liquids. 
Many theories rely upon this assumption. Therefore, in- 
formation about mobility of atoms is of great importance. 
In the case of glass-ceramics^, for example, one has to 
control nucleation and growth in the supercooled liquid 
to obtain the desired outcome. 

Despite of its importance, however, diffusion coeffi- 
cients are not always easily accessible through experi- 
ments. One largely employed method is to use tracers, 
which are atoms that can be monitored while they move 
within the material analyzed. The problem is that iso- 
topes that can be used as tracers are very expensive. A 
cheaper alternative is to employ an effective diffusion co- 
efficient obtained via the Stokes-Eintein or the Eyring 
relation. They state that diffusion divided by tempera- 
ture is inversely proportional to the viscosity of the liq- 
uid. Both relations differ only by a multiplicative factor. 
Determining viscosity is a much easier task and data for 
several materials are abundant. 

There is still an ongoing debate about the range of va- 
lidity of the Eyring relation. The relation is suitable for 
ordinary liquids above the melting temperature. As we 
approach the glass transition region, dynamical hetero- 
geneities in the liquid emerge and a novel diffusion regime 
takes placeS. Molecular dynamics have proven to be very 
useful to describe atomic processes on a nanometric scale. 
The following questions will be discussed in this work: 1- 
How can molecular dynamics assess the range of validity 
of the Eyring relation? 2- How do atoms rearrange in 
the system near the glass transition? 3- Are dynamical 
heterogeneities directly responsible for the breakdown of 
the Eyring relation?. 

According to recent studies, the breakdown of the 
Eyring relation depends mostly on the material's 
fragility. Studies on silica-, a strong glass former, re- 
veal that the Eyring relation holds all the way down to 



the calorimetric glass transition temperature Tg. On the 
other hand, it is shown for fragile liquids^ a breakdown 
near l.lTg. Simulations on model glass formers such as 
the binary Lennard-Jones liquid also point to a break- 
down. This breakdown also occurs far from Tg accord- 
ing to a recent study on supercooled water confined in 
nanopores^i. 

In this work, classical molecular dynamics simulation 
will be used to model lithium disilicate melt. It is an 
archetypal fragile oxide liquid for which transport data 
are numerous. Simulation results of diffusion and relax- 
ation will be compared to the experimental data avail- 
able. Also, a numerical experiment devised to track 
molecular rearrengements was performed. This work 
aims to shed light on the aforementioned questions and 
to give alternative interpretations about computer simu- 
lation results. 



II. COMPUTATIONAL DETAILS 

Classical molecular-dynamics (MD) simulation is the 
main computational method employed in this work. It is 
based on integrating the coupled newtonian equations of 
motion for a set of interacting particles. The particle in- 
teraction used to model lithium disilicate is the Bucking- 
ham potential, with parameters developed by Habasaki 
et alM. The potential has the form: 
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where r is the interparticle distance, i and j are the 
species, q is the effective charge for coulombic interactions 
and A, B and C are fitting parameters. This interatomic 
potential has proven to describe structural properties in 
the melt and the amorphous phases in good agreement 
with neutron scattering experiments. Further details on 
the potential can be found in Ref. |Tl|. 

All computer simulations were performed with the 
LAMMPS package, which is developed by S. Plimpton 
and co-workers ^?^. A time step of 1.0 fs was used in the 
entire work. Periodic boundary conditions were applied 
in all directions to simulate the bulk material. The long- 



range forces due to electrostatic interactions were prop- 
erly handled using the Particle-Particle Particle-Mesh 
(PPPM) techniquei^, with a real space cutoff of 10 A 
and energy accuracy of 1x10"'^. 

A system with 3888 atoms was setup initially in the 
orthorhombic structure for Liz 81203^1!^ at 300 K. To 
obtain liquid configurations, we heated the crystalline 
system up to 3000 K in a constant-pressure ensemble at 
zero pressure. Finally, the hot liquid was cooled down 
at a desired temperature and allowed to thermalize in a 
constant-volume ensemble. Thermostatting was done us- 
ing the Nose- Hoover chain method^^. During all stages, 
the simulation box was kept orthogonal in order to prop- 
erly simulate the liquid state. 

The liquid properties were calculated after adequate 
thermalization. The self-diffusion coefficient Di for the 
j-element was calculated via the Einstein relatior>ii 
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where \rf{t)) is the mean-square displacement at a time 
t. The structural relaxation time Tq was calculated from 
the intermediate scattering function^ 



(3) 



Fs{k,t) = {cos[k-ifit)-f{0))] 
and via fitting of 



Fsik,t) =expi 

T, 
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where k is the reciprocal vector of the first sharp peak 
of the static structure factor and f{t) is the position of 
a particle at time t. The structural relaxation time is 
proportional to the shear viscosity 77 via 



V{T) = GooTciT) 



(5) 



where Goo is the instantaneous shear modulus and T 
is the temperature. Goo is weakly dependent of the 
temperature^^. Hence, Tq and 77 have a very similar tem- 
perature dependence. 

The structural relaxation is a measure of how long the 
liquid structure rearranges globally due to viscous fiow. 
In the case of all alkali-silicates, those rearrangements 
are due to the mobility of the network-formers, i.e. Si04 
tetrahedra. For that reason, Ta was calculated for Si 
atoms only. 



III. RESULTS AND ANALYSES 

The primary goal of this work is to study the relation 
between diffusion and viscosity using the Eyring equation 



(see Sec lHI B]) and analyze its range of validity. To do 
this, we had to calculate D and Tq as a function of the 
temperature with the highest degree of undercooling as 
possible. A great care was taken in those calculations. 
D in the diffusive regime is properly calculated if, and 
only if, the mean-square displacement is linear with t for 
a prolonged period of time. Calculating I? on a sub- 
diffusive regime {\r'^{t)) oc f,^ <1) only gives us an 
upper limit for D. Also, to obtain a good fit for Tq, Fs{t) 
should be at least below l/e^. We have calculated D 
and Ta for Si atoms down to 1200 K, to which it was 
necessary a simulation run of 40 ns. 

The experimental value of the liquidus temperature for 
lithium disilicate is 1360 K. The orthorhombic-liquid co- 
existence temperature^^ calculated in this work is equal 
to 1450±30 K, in good agreement with experiment. 
Hence, MD simulation, in this case, could only achieve a 
small degree of undercooling to calculate relaxation and 
slowly diffusing atoms. Nevertheless, effects of under- 
cooling such as non-exponential relaxation were observed 
in our work. We have not detected any signs of crystal- 
lization during simulations, even at the highest degree of 
supercooling. 



A. Diffusion 

The temperature dependence of the self-diffusion for 
all three elements in lithium disilicate is plotted in Fig. 
[T] Along with MD data (filled symbols) are experimen- 
tal data (open symbols) available for comparison. Si dif- 
fusion was measured experimentally by electrochemical 
determination of cations interdiffusivity^*^ . Conductivity 
data at high temperatures^! were used to calculate Li 
diffusion via the Nernst-Einstein relation^^. An excel- 
lent agreement is found in Li diffusion between MD and 
experimental values. For the silicon diffusion, absolute 
values agree satisfactorily. The temperature dependence 
also shows excellent agreement with experiment. Over- 
all, Fig. [1] displays the typical behavior for an alkali 
silicate melt. The system presents two distinct diffusion 
mechanismsS^: a) slow Si and O atoms, which move as 
an Si04 unit (i.e. the network formers), and b) the light, 
fast moving Li atoms, which permeate the network and 
act as the network modifiers. It is believed that the for- 
mer mechanism is related to viscous flow over a wide 
temperature range accessible to experiments. This topic 
will be discussed on the next sections. 



B. Relaxation and the Eyring relation 

The intermediate scattering function was obtained and 
Ta was calculated from the fitting of Eq. S) This value 
of T„ in the equation [5] results in an effective viscos- 



ity rj 



relax 



which allows direct comparison with experi- 



mental viscosity data. Using Ta and the experimental 
viscosity value at 2200 K in Eq. [SI we obtained Goo 
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FIG. 1. Self-diffusion coefficients as a function of the recipro- 
cal temperature. Open symbols are experimental data for Li 
(Ref. [13) and Si (Ref. [2^. 



equal to 17.5 GPa. In Fig. [H 77'^'='°^ and a fit of the 
Vogel-Fulcher-Tammann function of experimental viscos- 
ity data3^ are plotted against temperature. One can ob- 
serve that both properties are incompatible except at the 
high temperature regime T >2000 K. Below that regime, 
the relaxation-defined viscosity is always larger than the 
actual data. 
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FIG. 2. Temperature dependence of relaxation-defined and 
experimental viscosities. 



The parameter /? was calculated via fitting of Eq|4] (in- 
set of Figl2|). If /3 = 1.0 the relaxation is exponential. 
It was shown that the system has non-exponential relax- 
ation for all temperatures analyzed. At 2400K (lOOOK 
above the melting temperature), the exponent /3 is close 
to 0.8 and decreases with further cooling. It is generally 
accepted that non-exponential relaxation is due to spa- 
tial heterogeneities found in the liquid upon supercool- 
ing. Note that in this silicate model, no undercooling 
was necessary to observe such behavior. This is an indi- 



cation that the non-Arrhenian behavior of the viscosity 
found in fragile liquids far above Tg is related to spatial 
heterogeneities. 

It is shown in Fig. [3] the temperature dependence of 
the characteristic structural length A calculated via the 
Eyriirg relation3^ 
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where A is a constant. Thus, the Eyring relation does not 
hold at any temperature using diffusion and relaxation 
calculated via MD. This observation had already been 
reported by Bordat et al^ for a binary Lennard- Jones 
model, and Saika-Voivod et al^ for a high-pressure frag- 
ile silica model. On the other hand, a simulation study 
on a metallic glass model^^ had shown that decoupling 
between diffusion and relaxation could only occur if the 
system were subject to deep undercooling. In our sce- 
nario, the degree of undercooling clearly does not play a 
role in this kind of breakdown. 

The interpretation of A comes from the original concept 
of Eq|B]as a hydrodynamic model and it is equal to 67rro, 
where rg is the apparent hydrodynamic radius. However, 
as the temperature decreases, particles have their mobil- 
ity restricted to cages of decreasing volume and diffusion 
becomes mediated by atomic jumps. The association of 
the characteristic length A and the size of the cage is 
still on debate and we will restrict the discussion of jump 
diffusion in the next section. 
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FIG. 3. Temperature dependence of the jump distance ob- 
tained using molecular-dynamics data. 



If Tj^'^^^^ is replaced by the experimental viscosity, the 
Eyriirg relation is valid all the way down to 1400 K. We 
applied a single A — 4.2 Ato the experimental data to 
allow comparison with diffusion calculated directly via 
MD (Fig. S]). Indeed, there are several experimental 
studiea^ i^^'^" showing that the Eyring relation is valid 
over a wide range of temperatures. Some studies on 



oxide glasses^i^ suggest a breakdown of the Eyring re- 
lation near the glass-transition temperature. Unfortu- 
nately, straightforward MD simulation is not capable to 
calculate diffusion at such low temperatures for silicate 
glasses. Still, it is worth mentioning that calculating dif- 
fusion coefficients at temperatures as low as 1200 K for 
this system can be problematic due to the prolonged re- 
laxations found in these conditions. Nevertheless, the 
trends observed in Fig. |4] are remarkable given the ex- 
treme conditions of the simulation. 
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FIG. 4. The calculated Si-diffusion coefficient and the Eyring 
relation as a function of temperature. A single jump distance 
A =4.2 A was sufficient to render the Eyring relation valid all 
the way down to 1200 K. Experimentally measured viscosity 
data were used in this plot. 



C. Description of jump mediated diffusion 

Another important question on the dynamics of fragile 
glass formers is about the way molecular rearrangements 
occur. Down to 1200K, atomic motion becomes rather 
sluggish and straightforward visual inspection is ineffi- 
cient. Based on the technique of A. Voter— i^ devised 
to identify infrequent events in many-body systems, we 
setup the following procedure: 1- Starting from a equili- 
brated configuration at 1200K, an MD simulation is car- 
ried out in the NVT ensemble; 2- A series of atomic con- 
figurations are recorded, being each sample separated by 
some fixed time interval; 3- These configurations are min- 
imized using a standard conjugate-gradient algorithm; 4- 
Si atoms displacements are calculated between two ad- 
jacent minimized configurations. The distribution of the 
displacements calculated for all minimized configurations 
are plotted in FigE] 

In agreement with several studies, this model for 
lithium disilicate also presented dynamical hetero- 
geneities, i.e., a subset of particles of the same species 
move farther than another subset at a given time pe- 
riod. FiglS] shows that particles that moved less than 



O 



100 



10 



^■^ 




lOps □ 






40ps • 
lOOps A 








□•^ 


, 


□ 






1 

1 


dHa 






H' 






E 






m % ^A 






^ • aA 



0.3 0.6 0.9 1.2 

displacement (A) 



1.5 



FIG. 5. Distribution of particle displacements between two 
adjacent configurations. Time intervals of 10 (squares), 40 
(circles) and lOOps (triangles) were used. 



0.8 A are well distributed in the samples. Between 0.9 
and 1.2 A we found a plateau- like distribution of roughly 
one particle. For larger displacements, distribution goes 
down to zero as expected. The suggested picture is that 
a small but important fraction of Si04 tetrahedra per- 
forms jumps that make up for the viscous flow observed 
from the diffusion coefficient. The remaining particles 
move, on average, a small fraction of the bond length. 
These sluggish rearrangements end up creating room for 
the jump mediated diffusion to occur. 

Additionally, FiglS] provides a concise analysis about 
caging and hopping processes. The well-distributed dis- 
placements can be associated with atomic vibrations oc- 
curing in spheric regions of radius close to 0.2 A . Note 
that for the two largest time intervals, the distributions 
are quite similar. On the other hand, the distribution 
height in the plateau changes considerably with the time 
interval applied. That reinforces the fact that diffusion 
is due to particles that perform jumps of sizes close to 
the Si-0 bond length. 



IV. CONCLUSIONS 

In this work, we tackled the problem of defining ef- 
fective diffusion coefficients from alternative dynamical 
properties. We applied molecular dynamics simulations 
in the entire work. The interaction potential to model 
our simulated lithium disilicate melt has been reported 
to reproduce quite well structural properties of both liq- 
uid and glassy states. This work reports that dynami- 
cal properties are also well reproduced for melts between 
900 and 2400K. Given the reliability of the potential, we 
were able to perform an analysis of the diffusion mecha- 
nisms on a realistic model for lithium disilicate. Results 
show that there is a decoupling between the structural 
relaxation time and the experimental viscosity. Also, the 



Eyring relation fails when viscosity is considered as a 
proxy of the strucutral relaxation. Direct comparison 
between the diffusion obtained in simulations and the ex- 
perimental viscosity revealed the Eyring relation is valid 
on a strict sense at least down to 1400K. Finally, an anal- 
ysis of the mobility of slow-diffusing Si-atoms showed two 
very distinct moving patterns. That observation con- 
firmed the emergence of dynamical heterogeneities on 
the system in the undercooled regime. Our conclusion 
is that the much discussed dynamical heterogeneities do 
not lead necessarily to a breakdown of the Eyring rela- 
tion. Rather, a decoupling between experimental viscos- 
ity and structural relaxation time obtained via MD takes 
place. Our study suggests that relaxation mechanisms 
are still not fully accounted for in MD simulations. For 



now, the Eyring relation continues to be a reliable mea- 
sure of kinetic processes in oxide melts except near the 
glass transition. 
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